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-^ ■ Abstract 

^ ' We use multiscale perturbation theory in conjunction with the inverse scattering transform to 

sij . study the interaction of a number of sohtons of the cubic nonhnear Schrodinger equation under 

^— N ' the influence of a smaU correction to the nonhnear potential. We assume that the sohtons are all 

\^ , moving with the same velocity at the initial instant; this maximizes the effect each soliton has on 

^^ ' the others as a consequence of the perturbation. Over the long time scales that we consider, the 

^^ , soliton amplitudes remain fixed, while their center of mass coordinates obey Newton's equations 

with a force law for which we present an integral formula. For the interaction of two solitons 
with a quintic perturbation term we present more details since symmetries — one related to 
^ ' the form of the perturbation and one related to the small number of particles involved — allow 

the problem to be reduced to a one-dimensional one with a single parameter, an effective mass, 
"t^ ' The main results include calculations of the binding energy and oscillation frequency of nearby 

O .. solitons in the stable case when the perturbation is an attractive correction to the potential and 

of the asymptotic "ejection" velocity in the unstable case. Numerical experiments illustrate the 
accuracy of the perturbative calculations and indicate their range of validity. 

rS 
H , 
. p. . 1 Introduction 

This paper is concerned with the asymptotic behavior of solutions of the initial-value problem for 
the perturbed nonlinear Schrodinger equation (NLSE) 

i5iV' + ^92^ + |V'|V+p[V',r]=0, (1) 

subject to the initial condition ^(x,0) = V'ol^;) for certain initial fields 'ipo{x), in the limit when the 
perturbation term p[^,^*] becomes formally small. The unperturbed problem, when p[^,^*] = 
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in (HI), is well-known to be solvable [||] by an inverse scattering transform, one consequence of 
which is the existence of finite energy soliton solutions that are dynamically stable and robust with 
respect to collisions. The unperturbed NLSE arises in two different physical situations in modern 
optics [||. Firstly, for high-speed telecommunications, (|l]) describes the propagation of light wave- 
packets along an optical fiber. In this interpretation, t is the spatial coordinate along the fiber 
and X is the retarded time variable for the signal; accordingly the solitons of the unperturbed 
problem (and usually also the solitary waves of the perturbed problem, when they exist) are called 
temporal solitons. The suggestion by Hasegawa and Tappert in 1973 1^, ^ that temporal solitons, 
being immune to dispersion, might serve as bits in a high-speed data stream has since generated 
a large body of work, much of which is comprehensively reviewed in |^, § . Secondly, for photonic 
switching devices, dTj) describes the stationary envelope of monochromatic light waves in a planar 
waveguide under the paraxial approximation. Here x and t are both spatial variables, with t being 
the propagation direction and x being the transverse direction; accordingly the solitons of the 
unperturbed problem are called spatial solitons. 

In both of these applications, the cubic term IV'PV' i^i @ models the Kerr effect in which the 
refractive index of the material depends linearly on the local intensity of light. For weakly nonlinear 
effects, when intensities are not too large, this effect is dominant in isotropic materials like glass. 
This fact, along with the integrability afforded by neglecting ^[■0, ^*], makes the unperturbed NLSE 
one of the most important models in modern optics. 

Of course, real materials can have a complicated dependence of refractive index on intensity, 
for which the Kerr effect is only an idealization. Modeling such phenomena requires introducing 
corrections to the coefficient \ip\'^ in the cubic term of the NLSE. The perturbative term pltpjip*] 
might also include corrections related to higher-order dispersion, the Raman effect, self-steepening of 
pulses, etc. In this paper we will consider only the influence of higher-order nonlinearity on solitons 
of the unperturbed NLSE. For spatial solitons in photorefractive media, such a perturbation can 
be the main factor influencing propagation. In particular, we take the perturbation in (|l|) in the 
form of a quintic term 

p[V,r]=<TeVlV, (2) 

where e > is a small parameter and o" = ±1. 

In view of the possibility of using solitons as bits in optical fibers or dynamically controllable 
switches in planar waveguides, it is of some interest to determine the effect of such a perturbation 
on the solitons of the unperturbed problem. If one considers an initial condition 'ipo{x) that is a 
"snapshot" of a simple soliton solution of the unperturbed problem, then there are many approaches 
available to study the perturbed evolution. Because the unperturbed soliton is stationary in some 
Galilean frame, the main effect ofp[i{j, ip*] will be an adiabatic adjustment of the soliton's amplitude 
and phase parameters. This fact, together with the simplicity of the form of the soliton solution, 
means that direct perturbative methods can be used to study their slow evolution. In particular, 
variational methods and multiscale methods applied directly to (|l|) often give valid results. These 
perturbative methods are dynamical in origin and capture effects on finite but long scales. Other 
methods can be used to answer infinite time questions concerning the persistence of solitary waves. 
In fact in the presence of quite general perturbations solitary waves continue to exist for arbitrary 
e 1]^, P and these can be expressed in closed form in some cases Q. 

The presence of more than one soliton complicates the analysis. If the solitons are isolated 
then the field may be approximated as a sum of solitons plus a small error term, and the adiabatic 
coupling among the solitons may be calculated by several methods. Note that if the solitons are 
moving with respect to each other then they will always be in isolation except possibly for a short 
time. An early analysis of this kind was carried out by Gordon [0], who studied the exact two-soliton 



solution of the unperturbed NLSE for equal velocities. When the solitons are well-separated, there 
is an effective force between them (even in the unperturbed NLSE) that varies sinusoidally with 
their phase difference. This phase difference grows linearly in t if the solitons differ in amplitude. 
The force is therefore zero on average [10| and one expects periodic motion. This is a physical 



explanation of the mathematical fact that the intensity \^p\'^ of the exact two-soliton solution for 
equal velocities is a periodic function of t. An extension of this argument to perturbed problems was 
given by Ankiewicz |1^], who obtained a simple description of soliton interactions with the use of 
complex averaged potentials. Again, the essential assumption is that the solitons are well-separated 
in X, so that the field may be approximated as a sum of solitons. If the solitons are close to each 
other, nonlinear interference effects cause the field to adopt a form very different from the linear 
superposition of individual solitons, and therefore a different approach is needed. Often, one turns 
to numerics to study the interactions of solitons in various media (see, for example, [^, ^, |l4|) 
without the restriction of the solitons being isolated. 

In the scattering transform domain, where the dynamics of the unperturbed NLSE are trivial, 
a state in which two solitons are close to each other in x has the same spectrum as a state in which 
they are far apart. This suggests that for studying the influence of perturbations on multisoliton 
bound states (that is, several solitons traveling with the same velocity, represented by a collection 
of eigenvalues of the Zakharov-Shabat equations with the same real part) it is best to carry out the 
analysis in the transform domain using soliton perturbation theory 0, |l5|, ^]. With pltpjXp*] ^ 0, 
the evolution of the scattering data is no longer trivial, and thus the scope of possible dynamics in 
near-integrable systems like (||) is much greater than in the unperturbed NLSE, including effects like 
repulsion, attraction, and energy exchange among bound or colliding solitons. Other techniques 
that have been used to study these effects include the judicious use of conserved quantities |2|, 
variational methods ||l^, 17|, "equivalent particle" approaches [18, 19], and of course, numerics. 



In this paper, we use soliton perturbation theory to study perturbations of the nonlinear po- 
tential in (|^), for initial conditions V'o(3^) that are snapshots of multisoliton bound states of the 
unperturbed NLSE. With respect to treating the solitons in isolation, this is a worst-case sce- 
nario since in the unperturbed NLSE a tightly-bound state of solitons will remain so for all time. 
Nonetheless, it is a scenario of some interest, in particular for the quintic perturbation @. If 
a = +1, then it is known that the solution remains bound, and this case has been studied using 
conservation laws [^]. If o" = —1, then the bound state becomes destabilized. Recently it was 
shown 1 21 1 by simulations of (|lj) that the instability causes the bound state to divide into isolated 



solitons that are ejected from the origin with nonzero relative velocities. On the time scales over 
which this splitting occurs, the solitons do not appear to exchange energy. In mathematical terms, 
each eigenvalue in the bound state ensemble, originally confined to the imaginary axis (zero veloc- 
ity), appears to slowly "grow" a real part while its imaginary part remains fixed. Once the solitons 
escape, they no longer interact and the velocities no longer change. The wave guidance properties 



of Y-junctions engineered from such splittings of spatial solitons have also been analyzed |22|. 

By considering the relative velocities to be small, we will find an integral formula that expresses 
the asymptotic velocity difference between a pair of initially co-propagating solitons destabilized by 
the quintic perturbation ^ with a = —1. Along the way, we will write down a coupled system of 
differential equations that describes the interaction of any number of solitons under more general 
perturbations over long time scales. These equations are just Newton's equations for a system of 
interacting particles in one space dimension; the particle coordinates have the interpretation of 
the soliton centers of mass. The force is translationally invariant, conserves the total momentum, 
and is also proportional to a, so the forces giving rise to attraction and repulsion are related just 
by a change of sign. For the interaction of two solitons, the problem may be reduced to a single 



degree of freedom, the relative separation of the solitons. The force law scales simply with the 
(fixed) amplitudes, which have the interpretation of masses. The result is a one-parameter family 
of problems indexed by a normalized effective mass. If the separation is small in the attractive 
case (T = +1, the force is nearly linear and the frequency of motion becomes a function of the 
normalized effective mass. We calculate this frequency, a quantity that is connected with the 
vibrations of solitons that are infinitely close, a limit opposite to the well-separated case. 

Our paper begins in §|| with a review of the theory of the scattering transform for the Zakharov- 
Shabat eigenvalue problem and of the inverse theory that holds in the reflectionless case. We also 
recall the derivation of the exact equations of motion in the transform domain corresponding to the 
perturbed NLSE (|l|). Then, in §^ we consider perturbations of the form p[i/j, tf;*] = e^W{\ip\'^)'ilj and 
apply multiscale perturbation theory to find asymptotic solutions of the equations of motion in the 
transform domain. The approximations are uniformly valid as e | on expanding time intervals 
of length e~^, and are given in terms of solutions of Newton's equations for particles interacting in 
one dimension under a force law that has several universal features. In §^ we focus on the quintic 
perturbation (^) and study the interaction of two solitons. We reduce the problem to the motion of 
a single particle and then explicitly perform the averaging required to remove secular terms from 
the asymptotic expansion. This leaves the force law in the form of a ID integral that we study 
numerically. We use it to compute the "ejection" velocity observed by Artigas et. al. [^] in the 
unstable case and the harmonic frequency of tightly-bound solitons in the stable case. Finally, we 
compare the results of perturbation theory with direct simulations of (|l|). The Appendix contains 
the more cumbersome formulae that nonetheless are among our main analytical results. 

Regarding notation, we will use stars for complex conjugation, and matrices will be written 
with bold letters, except for the Pauli matrices 
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2 Exact Inverse Scattering Theory For The Perturbed NLSE 

Here, we review the known inverse scattering theory for the Zakharov-Shabat eigenvalue problem 
to fix our notation. In general, we wish to consider (|l|) where p^ip, V'*] is & polynomial in ■0, •(/;*, and 
their x-derivatives. The field ip is taken to be in the Schwartz space as a function of x. 

2.1 Scattering data. 

We will work with the scattering transform of V'j a map that associates to the complex field ip 
at each fixed time a set of "scattering data" from which ifj can be reconstructed by inverting the 
map. As is well-known, the advantage of this is that the time evolution of the scattering data 
corresponding to the time evolution of ip is trivial when p = 0. Consequently, when \p\ -C 1, this 
proves to be a useful setting for perturbation theory. 

Fix f, and assume the complex function ■0(x,i) to be given. For A € M denote by M (a;,t, A) 
the 2x2 matrix solutions of the linear differential equation 



aM^ 



LM^ 



-iX -0 



M^ 



(4) 



satisfying the boundary conditions M (x,i,A)exp(iA(T3a;) ^ I as x ^ ±00 Since L is traceless, 
these boundary conditions guarantee that these matrices are unimodular for all x. For each A there 



can only be two linearly independent column vector solutions of (Q); therefore there is a matrix 
S(t, A), A € M, the scattering matrix, such that 

M'(x,t,A) =M+(x,t,A)S(t,A). (5) 

The first column of M~(x,t,A) and the second column of M^(x,f,A) turn out to be boundary 
values of analytic functions for S(A) > 0, while the second column of M~(x,t,A) and the first 
column of M^(x,t, A) are the boundary values of analytic functions for Q'(A) < 0. Adjoining the 
second column of M"'"(x, i, A) on the right of the first column of (|5|) and taking determinants gives 

Snit,X) = det{M- {x,t, \), M+ {x,t, X)) , (6) 

which is therefore the boundary value of a function analytic for ^(A) > 0. Likewise 5*22 (i, A) = 
det(M{'"(a;, t, A),M2~(x, i, A)) is the boundary value of a function analytic for Q{X) < 0. 

Fix A e M. Then, from (|), M=^* = a2M^a2, and thus S* = (J2S(T2, so that ^22 = S^^ and 521 = 
— 55^2- III particular, this means that as an analytic function for S(A) < 0, 522(i,A) = S'ii(t, A*)*. 
Also, for A G M the fact that det(S) = 1 implies the normalization condition ISnp + |5'i2p = 1- 

The analytic function Suit, A) with 9(A) > may have zeros Ai(t), . . . , AAr(t). The determinant 
formula @ then shows that there exist complex numbers 7i(t), . . . ,7Ar(t) such that 

M2+(x,i,Afc(t)) = 7fc(t)Mf(x,t,Afc(t)), k = l,...,N. (7) 

The conjugation symmetry of M^(x,t,A) for A G M, when extended to the complex plane, im- 
plies that at the complex conjugate points Xk{t)* where 5'22(i,A) vanishes, M^{x,t, Xk{t)*) = 
--fk{t)*M^{x,t,Xk{t)*), iov k = 1,...,N. Since 5ii(t, A) ^ 1 as A ^ 00 with 9(A) > 0, Hilbert 
transform theory can be used in conjunction with the normalization condition to express 511(4, A) 



for 9(A) > in terms of its zeros and the values of Si2{t, A) on the real axis |23]: 



N 



A-Afc(t) \ / 1 p log(l-|Si2(t,/i)| 



|2^ 



«"('. *> = (n jTi^) »p (^ /_ jzr, "^j ■ (8) 

The so-called "trace formulae" that equate certain functionals of the potential ip to functionals 
of the scattering data will be useful below. In particular, we will use the formula: 

/OO 1 /'OO _ _ 

%i,d^r)dx = — ^log(l-|5i2(t,;u)|2)dM-^9(Afc(t)2). (9) 

-00 ^"^ J —00 , _-, 

This functional (not to be confused with the perturbation p[V')V'*]) has the interpretation of the 
total momentum of the wave function ^(a;,t). For the unperturbed problem, as well as in the 
presence of many physically important perturbations, the total momentum is a constant of motion. 

2.2 Reconstruction of the potential in the reflectionless case. 

The miracle of inverse scattering theory is that for each fixed t, the potential tl^{x, t) can be recovered 
from its scattering data, namely the refiection coefficient Si2{t, A) for A € M, the eigenvalues {Afc(i)} 
with 9(Afc) > 0, and the proportionality constants {'^k{i)}- The reconstruction is particularly simple 
if 5*12 (i. A) = as a function of A for some i, since it then follows from ^ that 



which extends to 9(A) < as a meromorphic function. Similarly one sees that 522(i, A) = 
l/5ii(t, A) and that S'2i(i,A) = 0. Since S(t, A) is diagonal in this case, the solution matri- 
ces M^(a;,i,A) can be expressed in terms of a common solution matrix U(x,t,A) by setting 
M±(x,t,A) :=U(x,i,A)N±(t,A), where H 



N 



N 



l =Fl 
2 



N^it,X)=a,' diag [ll{X-Xk{t))-\l[{X-\kitrr']<T; 

\k=l k=l ) 

The columns of U(x, t, A) = (C/i(x, t, A), C/2(a^, ^, A)) necessarily satisfy the relations 

?72(x,t,Afc(t)) = 7fc(t)C/i(x,t,Afe(t)), -7fe(t)*C/2(x,f,Afe(t)*) = C/i(x,t,Afc(t)*) 

for all A; = 1, . . . , A^. It follows that U(x, t, A) takes the simple form 

/ ;v-i \ 

U(x,t,A) = A^I+ ^ APU(P)(x,t) exp(-iAa3x) , 



(11) 



(12) 



(13) 



that is, a polynomial in A times an exponential, where the matrix coefficients U''^''(x,t) are deter- 
mined uniquely from (|l2|). This means that (O) can be viewed as a linear algebraic system of \'N 
equations in 4A^ unknowns, the matrix elements of U^^^(x,i). Moreover, it can be shown that U 
constructed in this way is satisfies d^ = LU if and only if the potential function in L is 



i,{x,t) = 2iJJ^^-'-\x,t) 



(14) 



This formula reconstructs tlj{x,t) from the discrete scattering data {Afc(t)} and {'yk{t)} in the 
"reflectionless" case when 5i2(i, A) = 0. This treatment of multisoliton potentials via the matrix 



U follows Krichever [^], Manin |2g], and Date [p7| . See |28| for a relevant application. 



2.3 Dynamics of the scattering data. 

We now recall how the data evolve in t when if; satisfies (||). The motivating observation ||8| is that 
dl]) can be cast in matrix form: 



idth- d^B + [L,B] + P = 0, 
where the matrix L is the one appearing in the linear scattering problem 



(15) 



and where 



B 



A^ 



IV'P 



1 



-ixr-\d.r -A^ + ^iV'P 



plip.ip* 



(16) 



Using the fact that M satisfies (^, multiply ( p!5|) on the right by M and find 

{d^ - L){idt - B)M^ + PM^ = . 



(17) 



This equation is solved for {idt — B)M by variation of parameters. Introducing a new unknown 
J^(x,t, A) defined through the relation {idt — B)M^ = M^J^, one finds that J^ satisfies dxJ^ = 
— M PM . We now integrate to find J explicitly, taking into account the boundary conditions 



satisfied by M^ as x — > ±00 and the fact that in both hmits B — > X'^a^. With the use of these 
exphcit formulae for J the equations {idt — B)M = M J become equations of motion for the 
matrices M^: 



{idt - B)M=^ = M=^ ( -AV3 + / M=^"ipM=^ dx'] . 



(18) 



As written, (18) does not make sense for S(A) 7^ 0. But for A G M, the columns M^ and M2 
are the boundary values of functions analytic for ^(A) > 0, and we will also need equations for 
them that hold for 9(A) > 0. To this end, we introduce the matrix M(x,t, A) := {M^ ,M2), and 
as before define the new unknown J{x, t, A) = ( Ji, J2) through the relation (idt — B)M = MJ, and 
then integrate: 



J^ 



-A2 





M~ipMf dx' , J2 




A2 



+ 



M-^PM+ dx' . 



(19) 



As before, these expressions are used in {idt — B)M = MJ to yield the equation of motion for M, 
valid for 9(A) > except at {Ajt}, where M fails to be invertible. Each singularity is, however, 
removable, since detM = Sn and hence (writing M.^' for M.^(x',t,A)) 



M{x,t,\)M{x',t,\) 



-1 



Sn 



+1 



M7,M. 



MT'^M^' Mt,M- 



- A/r+' 



'12^"21 



12^" 11 



MfiMi2 



M2-1M2+2' - M2-/M2+2 Mf/M2+2 - M+2'^2"l 



(20) 



We make the natural assumption that the (isolated) zeros A = \k{t) of the denominator 5'ii(t, A) 
are simple [23|. But then the numerator of each entry is analytic at A = Afc(t) and is easily seen to 
vanish there, thus cancelling the singularity. Hence, the evolution equation for M makes sense as 
A — > Afc(i). We accordingly introduce the notation 



Hfc(x,x',t):= lim M(x, t, A)M(x', t, A)~^ 



(21) 



The equations of motion for M and M determine the evolution of the scattering data. Using 
S := M^ -"^M^, for real A one finds 

idtS = -M+-ifatM+ • M+"^M" + M.'^-^idtM.- = -M+-^i9iM+ • S + M+'^f^iM" . (22) 

Substituting from (ffl) yields 



idtS = AV3S - / M+ -1PM+ dx' -S- A^ScTs - S / M" "^PM" dx' . 

J X J —00 



(23) 



Finally, since S does not depend on x, it may be brought inside the integrals. With the use of its 
definition the integrals are combined, giving the equation of motion: 



/oo 
M+(x', t, A)-^PM- {x', t, A) dx' = 0. 
-00 



(24) 



Note that since P is off-diagonal, the equation for 5ii (t, A) only involves quantities analytic for 
9(A) > 0. Likewise, the equation for S22{t, A) only involves quantities analytic for 9(A) < 0. 
The equation of motion for the reflection coefficient 512 (t, A) is contained in that for S: 



/oo 
[M+(x', t, X)-^PM- {x', t, A)] ^2 dx' = 0. 
-00 



(25) 



The integrand here is p[^, V*]-^22-^22 ~P[^^ V'*]*-^]^-^]^) evaluated at x', t, and A, which generahy 
only makes sense for A € M, as required. Now, the expression defining the zeros Xk(t) of S'ii(t, A) 
is 5ii(t, Afc(t)) = 0. Differentiating with respect to t gives 

idtSii{t, \k{t)) + i'^it) ■ dxSnit, Xk{t)) = . (26) 

Using the equation of motion for S, one therefore finds 

The integrand here is |?[V',^*]M^M^ — p[^,^*]*M{2M{^, evaluated at x', t, and A = Afc(i). As 
remarked above, this makes sense with S(Afc(t)) > 0. It remains to find an equation for {7fc(i)}. 
Differentiating the defining relation M2{x,t,Xk{t)) = ^k{t)Mi {x,t,\k{t)) with respect to t and 
using the evolution equation for M taken in the limit A — ;• A^ {t) , yields the equation of motion 



z^(t)-2AS(t) 



/oo 
Hfc(x, X, t)PM{' {x, t, A) dx'+ 
'OO 

(28) 
i [dxM+ix,t,X)-jkit)dxM^{x,t,X)] ^(i). 



with A = Afc(t). The equations (|25D , (27), and (^Sj) describe the evolution of the scattering dat 



ta. 



but are coupled to the equations for M and M^. This coupling disappears for P = 0: 



idtSuit, A) - 2X^Su{t, A) = , i^{t) = , i^{t) - 2X,{tfjk{t) = , (29) 

for A: = 1, . . . , A^, as was first observed by Zakharov and Shabat |j|]. From this simple system, it is 
possible to introduce the coupling perturbatively, leading to closed systems order by order. 

3 Perturbation Theory with Nearly Bound SoUtons 

We now suppose that p[il>,ijj*] = e^VF(|^p)^ for some real-valued function W, taking e > to be 
a small parameter, and seek a perturbative solution of the equations of motion for the scattering 
data. We want a description of the solution up to an O(e^) error, containing important physical 
information, and valid uniformly over time scales of length 0(e^^). The initial data we consider is 

Si2(0,A)=0, Xk{0)=imk, 7fc(0) = exp(-2mfeX^ + i^^) . (30) 



Proposition 1 The solution of the initial-value problem of (p^, Wv^ ^'^^ (^^) with initial con- 
ditions (fgqj, is given asymptotically for small e by Si2{t, A) = O(e^) and 



Afc(i) = -^Vkiet) + imk + ©(e^) , 7fc(t) = exp(-2mfcXfc(et) + iOkit) + 0(6^)) , (31) 

where x^iT), v^iT), and 9k{t) are certain functions to he specified below. They satisfy Xfc(O) = x^, 
^fc(O) = and 0^(0) = 9^. This approximation is uniformly valid for times t = 0(e^^). 



We develop the expansion using the multiscale formahsni. Introducing the slow time variable 
T = et, and assuming all quantities to depend functionally on both t and T, we replace the time 



derivatives in (psj), (P7D, and (^) according to the chain rule: dt —> dt + edx- Observe that for 



the initial conditions (30), there is no enforced magnitude for K(Afc) or S'i2(A). We may thus select 
the scaling of these quantities to achieve a dominant balance. We choose to scale 5'i2(A) as e^ and 
5ft(Afc) as e. Thus, setting Afc = eak + ibk and 7^ = exp(Afc + i(,k)i we assume the expansions: 



euk + ibk = 


- e{af + ea^^^...)+^ibf+eb^H 


S12 = 


- e^(5g) + .5« + ...), 


^k + i^k = 


= (Am.A« + ...) + ^(er + ee« + 



(32) 



•) 



Substituting into the equations of motion and collecting powers of e, we examine the resulting 
equations order by order. First, from the leading-order terms in (^) we find for k = 1, . . . ,N that 

af' = ariT), bf=b^^\T), (33) 

so that these quantities do not depend on the fast time t. Similarly, looking at (^) we see that 

Af = Ai°)(T), er=er(0)-26r(r)2t. (34) 

The description we desire will follow upon determining the T dependence of these leading-order 
quantities. The 0(e) contribution in the equation for bk, the imaginary part of (p7|), is 

dtb'i^+9Tb^k^=0- (35) 

If this equation for bj^ is to be solvable in the class of bounded functions of t, then bj^ must be 
independent of T as well as t. With the T dependence of 6^, dropped, ( ^5[) can be solved by taking 
bj^ = 0. This yields the simplest part of the claimed result, that 9(Afc) is described uniformly 
for t = 0{e~^) by bk{t) = ruk + O(e^), where the mk are constants. Since b^, = nik, this also 
determines the leading-order behavior of ^^ from (p^). Setting 6^ := ^j^ (0), we define 9k{t) as 
follows 

ek{t):=c!'^=el-2mlt. (36) 

At 0(e), equation ( p8[ ) gives 

d,A^^^+dTAf=4afhf^=AmkafK (37) 

Again, avoid secular growth of A^ by setting 

5TAf(r)=4m4°)(r), (38) 

and then take A^ =0. If we now define: 

^kiT) := -^^ , MT) ■■= -2af{T) , (39) 



then (38) takes the simple form: 



An equation for Vk{T) is found at O(e^) in the real part of (|27|). We find 



dtal 



(1) 



-v'{T) = '^{Gk{t,T)), 



where Gk{t,T) is the leading term, divided by e^, of the right hand side of (^ 
from (^) and the leading-order behavior of {Xk}, we first see that 



dxSnit,Xk{t)) 



N 



*n 



A — irrii 



e=0 



X + irrij 
■1=1 J 



X=im]r 



1 T-r mfc - rrij 
2imk \}: ruk + rrij 



(40) 

(41) 
In more detail, 

(42) 



To find the leading-order behavior of M^, recall that S12 = O(e^) so that we can use the "reflec- 
tionless" construction of M and i/j in terms of U, which in turn is constructed from {Xk ~ iruk} 
and {7jfc ^ exp(-2mfeXfc(T) + i{^l ' (0) + 2mfct))}. This gives 

-1 



Gk{t,T) = i{-l) 



N 



2mfe JJ( 



^l 



m^ 



Wi\ijix,t)\'^)f{x,t,imk)dx, 



(43) 



with 



f{x,t,X) ■.= ij{x,t)U22{x,t,X)U2i{x,t,X) -ij{x,t)*Ui2{x,t,X)Uuix,t,X) . (44) 

Now, it is clear from (^) that all of the x and t dependence in U and ip enters through the 
products 7fc exp(— 2zAfca;) ss exp(CA;) exp(i6'fc(t)), where Ck '■= 2mfc(x — Xk{T)). Therefore, Gk{t,T) 
is a multiperiodic function of t for fixed T. The A^ — 1 frequencies are independent of T, since all 
of the T dependence enters through the functions Xk{T). Secular growth of a^ (t) is avoided by 
choosing f^(T') to cancel the mean value of this oscillatory function: 



mfe4(r) = FkixiiT), ..., XNiT)) := -2mk (9(Gfc(-, T))) 



(45) 



where angled brackets denote averaging over t with T fixed. The force functions Ff. depend para- 
metrically on the masses m^. Equations ( ^0[ ) and ( |45| ) imply Newton's equations for a system of 
interacting particles of mass mk and coordinate Xk : 



mkx'liT) = Fk{xi{T), . . . , xn{T)) 



(46) 



It is easy to see that Fk{xi + dx, X2 + dx, . . . , xjy + dx) = Fk{xi,X2, ■ ■ ■ , xjy) so that the forces 
only depend on the relative coordinates. There is also a symmetry for (^) coming from the 
conservation of momentum that holds exactly (and thus to all orders of expansion) in (|l|) with 
p['0,^*] = e^W {\'ip\'^)ijj . This symmetry follows from the trace formula (P) and shows that the total 
force on the system is zero: 



N 



N 



J2Fk{xi{T),...,XN{T)) = ^mkx'l{T) = 0. 



(47) 



k=l 



k=l 



The dynamical system ( [461) describes the evolution of the scattering data. Since the reflection 
coefficient vanishes to second order on the time scales of interest, solutions of ([46|) can be used to 



build, at each fixed t, the A-soliton potential as in §^j. This allows a direct comparison between 
numerics for (|l|) and the predictions of (^). 
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4 Two Particles 

Consider the case N = 2. The aforementioned symmetries imply that the system takes the form 

mi4(T) = Fiix^{T),X2{T)) = -\f{x2{T) - x,{T)) , 

m2x'i{T) = F2{xi{T),X2{T)) = \F{x2{T) - x^{T)) , ^^"^ 

for some function F. The relevant quantity is then the relative distance y{T) := X2{T) — xi(T), 
which has the simple-looking equation of motion 

my" = F{y) , (49) 

where the effective mass is defined by m := 2 (rn,^ + m^ ) 

4.1 Writing down the force function. 

We begin our study of the force functions by simplifying the integrand in ( ^3|) to isolate terms that 
are exact rc-derivatives and do not contribute. In this context, consider the squared eigenfunction 
system implied by (§). Let M be any solution of S^M = LM, and define the quadratic forms 

(P:=MuMu, X-=M2iM22, rj := M11M22 + M^Ahi . (50) 

Then, these quantities again satisfy a linear system of equations 

dx4> = -2iX4> + ijT] , 9a;X = 2iAx - V'*^ , 9^,7/ = -2V'V + 2V'X • (51) 

Using the quadratic forms associated with U, / as defined by (^) is seen to be an exact x-derivative: 

/ = ^d^V = \d^{UiiU22 + U12U21) = d^{Ui2U2i) , (52) 

where the last equality follows from the fact that the determinant of any solution of (Q) is indepen- 
dent of X because L is traceless. For A^ = 2, we use the relations ( ]l^ ) and the parameters Ai = imi, 
A2 = ir?i2, 71 = exp(— 2miXi(T) -|- i{9i — 2?7ift)), and 72 = exp(— 2?tt,2X2(T) -|- i(02 ~ 27712^)) to find 
U12 = e*^^(AV'/(2i) + if) and U21 = e-'^"" {\ip* / {2i) - ip*), where 



^ _ 2{ml-ml) 



D{Ci,C2,02-ei)) 



mi cosh(C2)e*^i(*^ - ma cosh(Ci)e*''2(t) 



V> 



(53) 

mim2(m|-mf) r . ^g u. . ^g ,^.i 

T^(. . n T\ smh(C2)e ^^ ' - smh(Ci)e ^^ > 

where ip is the well-known two-soliton "breather" solution, and using 

D{Ci,C2,(^) ■= {rrii -|-m2)^cosh(Ci)cosh(C2) - 2mim2 cosh(Ci -|- C2) - 2mim2 cos{9) . (54) 

Since W{-) £ M, only ^{f{x,t,imi:)) is needed to find 9(Gfc(t, T)). From (|5^ one finds / = 



-d^{X'^\ijj\'^ /A + X'^i'il^ip*) + lifl"^), and therefore 3?(/(x,t,imfc)) = mld^\ip\'^ / 4 - d^lipl"^ . Using this 
in the formula ( |4^ ) for 9(Gfc(t, T)), one finds that the first term is an exact derivative of a rapidly 
decreasing function and hence integrates away. In terms of the two quantities IV'P and |(/9p obtained 
directly from (53) we finally obtain 

9(Gi(t,r)) = — — 4 2T r Wm')dM^dx = -^Q{G2{t,T)). (55) 

2mi(m2 — m|) J_^ mi 
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In particular, it follows that —2miQ{Gi{t,T)) — 2m2''^{G2{t,T)) = so that the total instantaneous 
(that is, before averaging over t) force vanishes. 

Specializing further to the quintic perturbation (Q) by taking W{p) := ap^ and writing 



1 



27r foo 



F{y;mi,m2) = -T- / h{y,z,0;mi,m2) dzd9 , (56) 



2vrjo 
we have found the following explicit formula for h: 

, / /. N 128(Tm?m2(m2 - m?)^ „ , , , ^, 

h{y, z, 0; mi, ms) = ^/\ ^ [h + . . . + his] , (57) 

where the individual terms h^ are given in the Appendix. They depend on a dummy integration 
variable z that differs from x by a simple translation. Note that, by the periodicity with respect 
to the "fast" function 6 = 02{t) — 0i{t), averaging over t is equivalent to averaging over 0. 

4.2 Scale invariance. 

From ( |56[ ) and the explicit formulae for the terms h^ in the Appendix, note the important symmetry: 

F{^y; mi, 1712) = ^~^F{y; ^mi, ^7712) , (58) 

for all nonzero .^ € M. Setting y = Cq and S = ^~^T, the equation of motion takes the form: 

{^m)q"{S)=F{q{Sy,^mi,^m2). (59) 

For arbitrary masses mi and m,2, we may then set ^ = {mim2)^^'^. Because rh is homogeneous of 
degree one in mi and m2, it is convenient to use the normalized masses 

Ml = ^mi , M2 = ^m2 , M = ^m, ^ = (mim2)"^/^ . (60) 

Here, Mi and M2 satisfy M1M2 = 1 and may therefore be expressed in terms of the normalized 
effective mass M by solving 2{M2 + M^ )^^ = M subject to this constraint to find: 



Ml 



1-(1-M2)1/21 -M-^ M2 = [1 + (1 - M^)^/^l • M^S (61) 



assuming without loss of generality that M2 > Mi. From now on, we will work exclusively with 
the normalized masses, in which case the force depends only on S and M. 

4.3 Averaging. 

We now compute the ^-averages explicitly by residues. There are five terms: 

1 r^TT pQgp Q 2-7 r^^ cosP 

^Jo DiCi,C2,0y ~^Jo {a-cos0y 



\--=^l r.^7 / n^7 d0=^— I ,„ l:,,, dg, (62) 



for p = 0, 1, . . . , 4, where a := (2 — M^) • cosh(Ci) cosh(C2) ■ M~^ — sinh(Ci) sinh(C2) > 1- Changing 
variables to w = exp(i^), the contour of integration becomes the clockwise-oriented unit circle 
in the if-plane. The only singularity within the contour is a seventh-order pole at the point 
wq = a — (a^ — 1)^'^, where from here on the positive root is taken. Therefore, 



1 Res w^{w + w~^)P 

2P w=wo [w - woY {w - Wq^Y' 



A, = -- ""^ " ^",^ " >\^ . (63) 
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In particular one finds exact expressions for Ap = 65536(o^ — l)^^'^^p:, 

io = 8{2af + 240(2a)^ + 720(2a)2 + 160 , ii = 56{2af + 560(2a)3 + 560(2a) , 

A2 = 4(2a)6 + 232(2a)^ + 808(2a)2 + 192 , A3 = 42{2af + 588(2a)3 + 672(2a) , (64) 

A^ = 3(2a)'^ + 202(2a)^ + 928(2a)2 + 256. 

These results yield an explicit formula for the two-particle force function as an integral 

/oo 
Hiq,z)dz, (65) 

-00 

where we are assuming that M1M2 = 1 and M2 > 1 > Mi > 0, and where 

2o-fl - M^l^/^ ^ 
^ = £no — Yl [9mn{Mi,M2)tanhCi+g„m{M2,Mi)tanhC2]Hmn, 

m,n=0 

sech 2(6-"^) ^1 sech 2(6-") ^3 (66) 



m,n=IJ 



Hmn ■- . ^ ^ ^ 13/2 

-^ tanh Ci tanh (2 ) — sech ^(^i sech ^(^2 

Here, Ci = 2Mi(2: + q/2) and C2 = 2M2{z — q/2). Many of the coefficients gmn{oi,l3) vanish 
identically. In particular, gge = as is needed for the integral to converge. The nonvanishing 
coefficients gmnict,P) are given in the Appendix. 

4.4 General features of the force function. 



Unfortunately, ( pSf ) cannot be evaluated in closed form because the integrand generally involves 
both exp(Ci) and exp((^2)- Even if M2/M1 € Q so that the integrand becomes a rational function 
of, say, exp(Ci), the denominator is irreducible already for the simplest resonance, M2 = 2Mi. 
In spite of these difficulties, certain elementary features of the force law can be extracted: 

• F{q; M) is proportional to the constant a = ±1, as is clear from (|65|). 

• F{q;M) is an odd function of q, since the integrand satisfies H{—q,z) = —H{q,—z) and 
moreover this symmetry holds term by term in the formula for H. 

• F(q; M) decays to zero for large q. This follows from the fact that the denominator of each 
term Hmn in the integral is bounded and the corresponding numerator vanishes for large q 
whenever gmn 7^ 0. The result then follows from a dominated-convergence argument. 

• F{q; M) only vanishes exactly for g = 0. Thus it is strictly of one sign for g > 0. 

• The normalized effective mass M enters the dynamics both as a mass parameter multiplying 
the acceleration q"{S) and as a parameter in F{q;M) itself. 

The force F{q;M), as computed from the integral formula (|65|), is plotted in Figure |l| for several 
different values of the normalized effective mass M. 
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Attractive Case 




-5 

Displacement q 



Figure 1: The force law F{q,M) in the attractive case, a = +1, for M 
repulsive case a = —1, the force simply has the opposite sign. 



0.3,0.4,0.5,0.6. In the 



4.5 Attractive case. Spring constant. 

For a = +1, the force F{q;M) and the displacement q have opposite signs, so the force is always 
attractive. This means that the slow dynamics of the two-soliton bound state are periodic in time 
and the state remains boundg. To illustrate. Figure y compares the results of perturbation theory 
with a simulation of ([l|). For small displacements, we have F{q;M) = —k{M)q + 0{q^). The 
(mass-dependent) spring constant k{M) determines the frequency uj{M) := {k{M)/ MY'"^ of small 
oscillatory motions. This is the frequency on the time scale S; the frequency on the original time 
scale t is related by r2(r?2i, rn-2, e) = e {mim2)^''^uj{m/{mim2Y''^). A formula for the spring constant 
k{M) can be found by simply differentiating with respect to q in ( |65|) and setting (7 = 0, however it 
seems less useful to present than a plot, shown in Figure y, of the (numerically) evaluated formula. 
In Figure |^ we plot the corresponding frequency (on the time scale 5), the latter being a directly 
observable quantity. It is noteworthy here that the dynamics of solitons can be described by a 
linear theory even though their amplitudes are not at all small. The parameter linearizing the 
theory is the distance between the solitons, rather than the soliton amplitude. We also remark that 
the limit in which this linear behavior holds is that of infinitessimally-separated solitons, a limit in 
which methods assuming the solitons to be well-separated are invalid. 



4.6 Repulsive case. Asymptotic velocity. 

For (T = — 1, the force and displacement q have the same sign, resulting in q always becoming large. 
Solitons that are near each other at t = are ejected from the origin as observed by Artigas et. 
al. |21]. This effect is captured accurately by our theory, as shown in Figure ||. The work done by 
the force in moving the particle from q = qQ to q = oo determines the asymptotic velocity of an 
initially stationary particle upon ejection. Taking q^ = {) corresponds to the ultimate velocity of a 



^This is a long-time statement, holding for t — 0(e ^), but not an infinite time statement, 
whether true breather-like bound states exist (that is, permanently) for nonzero e is more subtle. 



The question of 
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Figure 2: A two-soliton bound state affected by an attractive perturbation. Here, e = 0.0387, 
mi = 0.6, and ni2 = 1. Left: an approximation to IV'P found by solving Newton's equations and 
then constructing the field using reflectionless inverse theory. Right: the corresponding numerics 
for fUj. The bound state has too much energy for the harmonic approximation to hold, and the 
period of motion, about 120 time units, is longer than the harmonic period. 



stationary particle that is slightly perturbed from (unstable) equilibrium at the origin. With zero 
initial velocity, one equates the asymptotic kinetic energy with the work done: 



\Mq'{^f 



F{q-M)dq, 



<?o 



to find a formula for the asymptotic velocity difference: 



q'{oo) = (i- 



H{q, z) dz dq 



1/2 



(67) 



(68) 



Figure ^ shows the asymptotic velocity difference q'{oo) for go = found from (|^) as a function of 
the normalized effective mass M. To apply the graph in Figure ^ to problems with unnormalized 
masses, it is useful to unravel the changes of variables made so far. Given mi and m2, the scaling 
parameter is ^ = (mi?Ti2)~^'^ and the effective mass is m = 2-(m^ +^^2" )~^ ■ Then, the normalized 
effective mass used in Figure |6| is M = ^m. Next, from the graph one finds the asymptotic velocity 
g'(oo). The true velocity in the original coordinates is then dy/dt = e^~^g'(oo). For example, the 
parameters used in Figure ^ imply a normalized effective mass of M ~ 0.9. From Figure ^ one 
finds g'(oo) ~ 5.0, and thus dy/dt ~ 0.15. This value agrees well with the pictures in Figure 5[ 

In the attractive case, the integral ( p7| ) also has physical meaning as the binding energy of the 
two-soliton state. A relative velocity in excess of g'(oo), the escape velocity, will "ionize" the state. 



5 Discussion 

Multiscale asymptotics shows that under certain conditions the behavior of a multisoliton initial 
condition in a perturbed NLSE reduces to Newton's equations for a system of interacting particles, 
one particle per soliton. The theory applies over time scales of length 0(e~^) for perturbations 
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Figure 3: The spring constant for small hound motions as a function M. 

of size e^, when the initial velocities of the solitons mutually differ by an 0(e) amount. Our 
calculations make very concrete the often-cited analogy between solitons and particles. We want 
to emphasize that the limit considered here is one in which the relative velocities of the solitons 
are small but the solitons may be strongly nonlinearly superimposed, precisely the limit in which 
methods exploiting large distances between solitons fail. 

For a quintic perturbation of the NLSE and an initial condition composed of two solitons, 
the resulting dynamical system can be analyzed. When the perturbation is attractive (cr = +1), 
the system describes a nonlinear oscillator with all solutions q{S) being periodic. If the energy 
associated with q{S) is small (that is, if q{0) and q'{0) are both small), then the periodic motion 
is nearly harmonic, and formulae for the associated spring constant and frequency of motion can 
be found; in this limit the model for the soliton interaction linearizes even though the soliton 
amplitudes are not at all small. The latter are determined by the masses mi and ni2 and are not 
related to the coordinate q{S). For larger energies, the spring "softens" and the frequency decreases 
with increasing energy. The pictures in Figure y show oscillations in the nonlinear regime, where 
the frequency of motion is smaller than the linear frequency. Of course even in the nonlinear 
regime, the dynamics still obey the simple model Mq" = F{q;M). Although the periodic motion 
is predicted and observed over long time scales of size 0(e~^), it is not likely to persist for all time, 
due to the influence of higher-order resonant coupling effects. 

On the other hand, when the quintic perturbation is repulsive {a = —1), the nodal point at 
the origin in the phase plane gets replaced with an unstable saddle point. All orbits apart from 
the fixed point itself represent the nonlinear development of the instability. Because the force 
vanishes fast enough for large q, the velocity q'{S) ultimately saturates as the two-soliton state 
becomes "ionized". From the force function F{q;M) this "ejection" velocity may be calculated, 
giving excellent agreement with direct simulations of the perturbed NLSE. This analysis explains 
the observations reported in [^ . The symmetry-breaking that determines which soliton ends up on 
the right and which on the left can be traced to the location of the initial phase point in relation to 
the separatrix connected to the saddle. Unlike in the attractive {a = +1) case, the approximation 
obtained from multiscale asymptotics for the repulsive (cr = —1) case is expected to be uniformly 
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Figure 4: The frequency lo of harmonic motion as a function M . 

valid for all time, since as the solitons separate, further effects due to resonant coupling diminish. 
Given the formula for the force F{q;M), it is possible to compute the harmonic frequency 
and ejection velocity, more explicitly than we have done here. For example, the formulae would 
be expected to simplify in the limits M J, (corresponding to two solitons differing very much 
in amplitude) and M j 1 (corresponding to two solitons with nearly the same amplitude). The 
calculation of the ejection velocity is challenging because it may require uniform approximation 
of F{q; M) for all q in the limit of interest; pointwise asymptotics for fixed q are not enough to 
approximate the work integral without further information. 
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Appendix: Formulae for the Two-Particle Force Function Integrand 

Here, we record the details of the formulae for the two-particle force function needed to calculate 
or approximate for special values of M the force and related quantities to any desired accuracy. 



Before averaging. The thirteen terms appearing in the sum in ( |57| ) are given here in terms of 
c := cos 6*, Sk '■= sinh(("fc) and Ck := cosh(Cfe), Ci '■= 2mi{z + y/2) and (2 '■= 2m2{z — y/2). 

hi = 2mlm2S2CiCl + 2mimlSiCfC2 /i2 = - (ml{ml + ml)SiCl + ml{ml + ml)S2Cl 



/i3 = 2mf (m? + mlc^)SiCl + 2m^(m?c^ + m^)52Cf 
/14 = -{2m\{ml + m\)cS2Cl + 2ml{ml + rn\)cSiCl 
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Figure 5: A two-soliton bound state affected by a repulsive perturbation. Here, e = 0.07746, mi = 
0.4, and m2 = 1. As in Figure]^ the result of perturbation theory is on the left and the numerics 
are on the right. The solitons escape with a relative "ejection" velocity given by (|g^. 



/is = mliml - 9ml)cS2CfCl + m^(m^ - 9m?)c5iC7f C| 

he = m\m2{^rn\ + 3m^)c5iCiC| + ?TT,i?n,|(3?n,2 + f)m\)cS2CfC2 

h7 = - (2mfm2{{5ml + m^)c + 4m^c^)SiCiC| + 2mim^((m? + 5m|)c + Amlc^)S2CfC2) 

hs = 2m\m2(m\ + (5m? + ^nil)'?)S2CxC^ + 2mim\{ni\ + {hm\ + ^m\)(?)SxC\C2 

hg = -(imlmiil + 4c^)SiCfCl + 2miml{l + Ac^)S2CfCi^ 

hio = mfm2{3ml - mj){l + 4c^)S2CfCl + mimiiSmj - m|)(l + 4c^)SiCfCl 



hu = -[2mfmi{mi - m|)(3c + 2c^)S2QC| + 2mfm^(mf - m^)(3c + 2c^)5iCfC7| 

hi2 = -Urnlrnl{{3ml+ml)c+{2ml+Aml)c^)SiC'fC^+Arnlml{{3ml+ml)c+{2ml+Aml)c^)S2CfCl 
his = 4mlml{ml + {3ml + 4ml)c^ + 2mlc^)S2CfC^ + Amlml{ml + {3ml + 4ml)c'^ + 2mlc'^)SiCfCl 

After normalization and averaging. Here, we give the nonzero quantities gmn = gmnict,l3) 
appearing in (|6q). In these expressions /? and a are Unked by the normalization condition a/3 = 1. 



503 
505 



672a3 - 672a'^ 
-2304a3 _ 2816a^ 



504 = 1344a3 + 3136a^ 
506 = 512a3 + 512a'^ 



512 
513 
514 
515 
516 



-60480^ - 40320^^ + 10080/3 
69664a3 + 16576a'^ + 11200^1 + 2240/3 
-130544q3 - 52016a'^ - 29600^1 - 29520/? 
759040^ + 48256a'^ + 44800^^ + 18816/? 
-10368^3 - 10368a^ - 1920a" - 1920/3 
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Figure 6: The asymptotic velocity difference q'{oo) of two solitons falling from unstable equilibrium. 



921 
922 
923 
924 
925 
926 



930 
931 
932 
933 
934 
935 

936 

940 
941 
942 
943 
944 

945 



= -13440q3 + 3360/3 + 10080/3^ 

= 17920a3 + 17920a'^ + 51520/3 - 24640/3^ 

= -2347200^ - 48320a^ - 54400^^ - 272640/3 - 7200/3^ 

= 479872q3 + 146176a^ + UlSAa^^ + 3200^5 + 352704/3 + 36864/3^ 

= -2771520^ - 144096q'^ - 19280a" - 560a^^ - 141808/3 - 15120/3^ 

= 31680a3 + 31680a'^ + 88000^^ + 480a^5 + 8800/3 + 480/3^ 

-4032/3 + 3360/3^ + 672/3^ 

5600a3 - 90944/3 - 32480/3^ - 7616/3^ 

1316000^ - 67200"^ + 166960/3 + 14960/3^ + 15760/3^ 



-24576q3 - 
-375744a3 



15168a^ + 3648a" + 249088/3 + 145344/3^ - 1984/3^ 
63840a'^ - 6688a" - 480a^^ - 545056/3 - 186144/3 



2803680^ + 1087680^^ + 13712a" + SUa"^ + IQa^'^ 

+ 227744/3 + 54604/3^ + 2592/3^ 
-24024a3 _ 24024a'^ - 8008a" - 728a^^ - 8a^^ - 8008/3 - 728/3^ - 8/3 

-6720/3 - 22400/3^ - 2240/3' 
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31680q3 + 216128/3 + 162400/3^ + 19072/3^ + 800/3 
-269952a3 - 187200^^ - 592576/3 - 268928/3^ - 29984/3^ - 2560/3 
384912a3 + 82032a'^ + 2976a" + 387040/3 + 43584/3^ - 8784/3^ + 1168/3 
-81312a3 - 63840a'^ - 8256q" - ima^^ 

+ 97152/3 + 127104/3^ + 26784/3^ + 864/3^3 
-61776a3 + 4368a" + 288a^5 - 96096/3 - 39312/3^ - 4032/3^ - 48/3 



13 



■13 
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550 = 9216;3 + 21760/3^ + 4864/3^ 

551 = -20096a3 - 129152/3 - 155264/3^ - 38272/3^ - 1280/3^3 

552 = 112848a3 + 7584a'^ + 355920/3 + 307872/3^ + 71872/3^ + 3984/3^3 + 80/3^^ 

553 = -165584a3 - 22752a^ - b2Sa^^ 

- 330528/3 - 2001 12/3^ - 32672/3^ - 1392/3^3 - 96/3^^ 

554 = 72072a3 + 15288a'^ + 648a^^ 

+ 92664/3 + 24024/3^ - 6552/3^ - 1512/3^3 - 24/3^^ 

560 = -1536/3 - 4096/3^ - 1536/3^ 

561 = 1280a3 + 13568/3 + 27648/3^ + 13568/3^ + 1280/3^3 

c,g2 = -4096a3 _ gea"^ - 27808/3 - 50688/3^ - 27808/3^ - 4096/3^3 _ gg^i? 
563 = 2912^3 + 112a'^ + 16016/3 + 27456/3^ + 16016/?^ + 2912/3^3 _^ 112/317 
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